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new computational framework for the simulation of turbulent flow through complex objects and 
^yi^^long irregular boundaries is presented. This is motivated by the application of metal foams in com- 
^^^act heat-transfer devices, or as catalyst substrates in process-engineering. The flow-consequences of 
• i-^uch complicated objects are incorporated by adding explicit multiscale forcing to the Navier-Stokes 
'-Equations. The forcing represents the simultaneous agitation of a wide spectrum of length-scales 
^^hen flow passes through the complex object. Two types of forcing procedures are investigated; 
^S^ffith reference to the collection of forced modes these procedures are classified as 'constant-energy' 
'constant-energy-input-rate'. It is found that a considerable modulation of the traditional energy 
I C ascading can be introduced with a specific forcing strategy. In spectral space, forcing yields 
strongly localized deviations from the common Kolmogorov scaling law, directly associated with 
I the explicitly forced scales. In addition, the accumulated effect of forcing induces a significant 
'^non-local alteration of the kinetic energy including the spectrum for the large scales. Consequently, 
^j-a manipulation of turbulent flow can be achieved over an extended range, well beyond the directly 
\^^^OTced scales. Compared to flow forced in the large scales only, the energy in broad-band forced 
^Urbulence is found to be transferred more effectively to smaller scales. The turbulent mixing of 
p, ' passive scalar field is also investigated, in order to quantify the physical-space modifications of 
■ transport processes in multiscale forced turbulence. The surface-area and wrinkling of level-sets of 
^~-^he scalar field are monitored as measures of the influence of explicit forcing on the local and global 
ixing efficiency. At small Schmidt numbers, the values of surface-area are mainly governed by 
Othe large scale sweeping-effect of the flow while the wrinkling is influenced mainly by the agitation 
""^af the smaller scales. 

. ^H^'^ paper is associated with the focus-issue Multi-scale Interactions in Turbulent Flows. 

^ : 

01 |. Introduction 

• '"jVarious multiscale phenomena in turbulent flows arise from the passage of 
r_jfl,uid through and along geometrically complex objects placed inside the flow- 
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domain. The corresponding perturbations of the flow arise simultaneously on 
a range of length-scales and find their origin in the complexity of the bound- 
aries of these objects. A motivating example is the flow through a porous re- 
gion such as a metal foam depicted in Fig. ^ Many more examples can readily 
be mentioned, arising in different technological applications or in numerous 
natural flows, including flow over forest canopies [1,3,12]. 




Figure 1. A porous nickel foam contains various geometrical complexities on different 

length-scales [29]. 

The purpose of this paper is to investigate the computational modeling of 
flows through complex regions via the introduction of explicit forcing terms 
in the Navier-Stokes equations. Consistent with the many shape-details of 
the obstructing objects, such forcing will need to represent the perturbation 
of the flow on various length-scales simultaneously. This distinguishes the pro- 
posed computational modeling from more conventional forced turbulence pro- 
cedures. In the latter the flow agitation is restricted to a few large scales only 
with the aim to observe the development of a natural inertial range at smaller 
scales in the turbulent flow [24,46]. Instead, in this paper we allow the forcing 
of a collection of widely different modes. The consequences for transport and 
dispersion in such turbulent flows will be studied both in spectral - as well as 
in physical space. We will primarily establish the degree by which the spec- 
tral properties of a turbulent flow can be modified relative to the classical 
Kolmogorov scaling, and quantify the efficiency with which embedded scalar 
fields can be mixed by the modulated flow. 

Complementary to the proposed explicit forcing approach, two alternative 
formulations have been put forward in literature to capture the flow in and 
around complex objects. These include the explicit boundary modeling [49] 
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as well as an approximation in terms of effective boundary conditions and 
(surface) roughness parameters [20,22]. The roughness parameterization has 
been introduced for situations in which the roughness length-scales are much 
smaller than the boundary layer thickness [44]. For geometries which display 
both large- and small-scale contortions of the shape of the object, compared 
to the boundary layer thickness, the surface-roughness parametrization may 
not be sufficiently accurate [10]. Alternatively, in case of explicit boundary 
modeling, no-slip conditions are imposed at all the intricate shape-details of 
the object. This computational approach can in principle achieve full accuracy 
but is limited to cases of modest complexity in view of the elaborate geometric 
modeling and the high computational expenses that are required [4,5]. 

Flow through complex gasket structures may give rise to self-similar turbu- 
lence spectra which do not follow the well known Kolmogorov —5/3 slope [27]. 
Such non-Kolmogorov turbulence was observed in flows over tree canopies, 
and is reminiscent of a spectral short-cut feature that was also observed exper- 
imentally [12]. In this paper we investigate the potential of multiscale forcing 
to accurately characterize such dynamic flow-consequences of complex domain 
boundaries without the need to explicitly account for their intricate geometri- 
cal shape. We consider the incompressible Navier-Stokes equations with mul- 
tiscale forcing working as a stirrer whose dynamical effects are controlled by 
a distribution of simultaneously perturbed length-scales. 

To arrive at a multiscale modeling that is quantitatively linked to actual 
complex objects several steps need to be taken. In this paper we address 
a flrst step in which we examine in some detail the influence different forc- 
ing procedures have on the energy dynamics in spectral space and the mixing 
characteristics in physical space. Special attention is devoted to the mixing 
efficiency of a passive tracer by monitoring the surface area and wrinkling of 
level-sets of these scalar fields [14]. Specifically we look at the instantaneous 
and accumulated effect on surface area and wrinkling caused by broad-band 
forcing. 

Different divergence-free forcing procedures will be applied to directly per- 
turb a large number of flow-scales. The alterations of the flow dynamics express 
themselves clearly in the kinetic energy. The transfer of energy toward smaller 
scales is found to increase considerably, compared to the case in which only 
large scales are forced. When a specific narrow band of scales is agitated by 
the forcing, then the locally higher spectral energy is not 'compatible' with 
the molecular dissipation-rate and an accelerated transfer is observed toward 
smaller scales. This effect is found for both families of forcing methods, i.e., 
constant-energy and constant-energy-input-rate. The kinetic energy spectrum 
is also modified non-locally, in a range of scales that are larger than the di- 
rectly forced scales. Consequently the agitation of a band of small length-scale 
features can accumulate and also induce significant alterations of the largest 
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flow- features, e.g., by contributing to an increased backscatter. 

Tfie clianges in the flow-dynamics due to the application of broad-band forc- 
ing application also has consequences for the turbulent transport properties 
of the flow. This may be expressed in terms of the mixing efficiency of embed- 
ded passive scalars. In particular, monitoring the surface-area of level-sets of 
the passive scalar allows to characterize changes in the large-scale sweeping of 
the flow, due to the forcing. Likewise, the more localized motions directly affect 
the 'wrinkling' of the passive scalar level-sets. The dependence of these mea- 
sures for the mixing efficiency on forcing parameters can be used to quantify 
the mixing efficiency arising from agitation of different bands of flow-structures 
with different forcing strengths. Specifically, we investigate the dispersion of 
strongly localized initial scalar concentrations. The direct numerical simulation 
of the forced turbulence shows that the maximal surface-area and wrinkling 
as well as the time at which such a maximum is achieved can be controlled 
by variation of forcing parameters. The time- integrated surface-area and wrin- 
kling are indicators of the accumulated effect. The simulations show that at 
small Schmidt numbers, a higher emphasis on small-scale flow agitation yields 
a significant increase in the time-integrated total mixing of the flow. 

The organization of this paper is as follows. In section [21 the simulation 
method, together with the explicit forcing strategies are introduced. Section |31 
is devoted to the modulation of the cascading process associated with the dif- 
ferent forcing methods. The consequences of forced turbulence for transport 
and dispersion in physical space will be quantified in section |1J Concluding 
remarks are collected in section |21 



2 Simulation of forced turbulence 

In this section we will first introduce the governing equations (subsection 12.111 
and subsequently describe the explicit forcing strategies that are used to drive 
the flow (subsection 12. 2|) . Two types of deterministic forcing strategies will be 
included: procedures which yield constant-energy in the collection of forced 
modes, and procedures which correspond to a constant-energy-input-rate for 
these modes. Subsequently, the computational method, its validation and par- 
allel performance will be described (subsection I2.3|l . 
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2.1 Governing equations 

The dimensionless system of nonlinear partial differential equations which gov- 
erns the flow of a viscous incompressible fluid is given by: 

^^^^ + (v(x,t)-v)v(x,t) = -^Vp(x,t) + i.v2v(x,t)+f(x,t) 
^ V-v(x,t) = 

where v is the velocity field, p is the constant mass-density and p the pres- 
sure. The dimensionless viscosity is the inverse of the computational Reynolds 
number Re, i.e., ly = 1/ Re, and f is the external forcing which we will specify 
in subsection 12.21 This system of equations may be rewritten in terms of the 
vorticity w(x, t) = V x v(x, t). Making use of the identity: 

(v(x,t) • V)v(x,t) =L^(x,t) X v(x,f) + iv(|v(x,t)p) (2) 

we may express as: 
d 



i/V^J v(x,t) =w(x,t)-V(^-p(x,t) + -|v(x,t)|^J +f(x,t) (3) 

where we introduced the nonlinear term w(x, t) = v(x, t) x u;(x, t). 

The flow domain is assumed to be periodic with the same period in each of 
the three coordinate directions. An efficient representation of the solution in 
terms of Fourier modes can be adopted [6,33,52] in which the velocity v(x, t) 
is expanded as: 

v(x,t)=^u(k,t)e^'^- (4) 



and the wavevector k (fc = |k|) has components ka = 2TTna/Lh, 
Ua = 0, ±1, ±2, . . . for a = 1,2,3. The dimensionless length of the periodic 
domain is denoted by Lf, and iiQ,(k, i) is the Fourier-coefficient correspond- 
ing to the k-th mode of Vq,(x, t). The equation governing the evolution of 
the Fourier-coefficients is given by: 

+ lykA u(k, t) = W(k, t)-ikT (-p{x, t) + ^ |v(x, t)\^,-k\ +F(k, t) (5) 



where .F(a(x, f),k) denotes the Fourier-coefficient of the function a(x, t) cor- 
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responding to wavevector k: 

^(a(x,t),k) = A(k,t) if a{x,t) = ^A{k,t)e'^'' (6) 



In addition, W(k, t) and F(k, t) denote the k-th Fourier-coefficient of the non- 
linearity w(x, t) and forcing f(x, t), respectively. 

In spectral space the pressure term may be eliminated from ((S)) if use is 
made of the incompressibility condition. This is equivalent to the well-known 
practice of solving a Poisson equation for the pressure in physical space for- 
mulations [48] . If we multiply © by k, use the continuity equation in spectral 
space, i.e., k • u(k, t) = 0, and assume that the forcing itself is divergence- free, 
so that k • F(k, t) = 0, the pressure term can be written as: 

-.-/I / N 1 , / m2 , \ k- W(k,t) 

.F(-p(x,t) + -|v(x,t)|2,kj = (7) 

The equation for the Fourier-coefficients of the velocity field © may now be 
written as: 

I + uk'^ u(k, t) = W(k, t) - k(^i-^^) + F(k, t) (8) 

This may be expressed in a more compact form in terms of the projection 
operator D defined by: 

D^, = 5^,-^ (9) 

This operator restricts the solution to the space of divergence-free fields, 
represented by Fourier-coefficients u(k, t) which lie in the plane normal to 
the wavevector k. We obtain the governing equation for the desired Fourier- 
coefficients as: 

^ + u(k,t) = DW(k,t) + F(k,t) (10) 

A more detailed discussion of this spectral approach to the Navier-Stokes 
equations is available in [33]. It forms the basis for the numerical treatment 
that will be specified in subsection 12.31 

In various applications the dispersion of a passive scalar by a turbulent flow 
is of central importance. Passive scalar transport may be used to characterize 
the physical space consequences of multiscale forced turbulence. The governing 
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equation for the evolution of the scalar concentration C(x, t) contains advec- 
tion by the velocity field v(x, t) as well as diffusion. In physical space this may 
be expressed as: 

+ (v(x,t) . V)c(x,t) = KV2c(x,t) (11) 

where k is the non-dimensional molecular diffusivity of the scalar. Compared 
to the dimensionless viscosity in we adopt k = v/Sc where the Schmidt 
number Sc characterizes the scalar diffusion. Roughly speaking, if Sc > 1 then 
the scalar field displays a wider range of dynamically important length-scales, 
compared to the turbulent velocity field, while values Sc < 1 indicate a com- 
parably smoother scalar field. The equation which governs the development of 
the Fourier-coefficients c(k, t) of the scalar field C(x, t) can readily be found as: 



(^■^ + Kk'^^ c(k, t) = Z{k, t) where Z(k, t) = ^((v(x, t) ■ V)C(x, t), k) 

(12) 

The changes in the turbulent transport properties of the fiow due to the mul- 
tiscale forcing can be investigated by considering the evolution of the scalar 
concentration at different Schmidt numbers. The structure of the left-hand 
side of (|12j) is identical to the Navier-Stokes equations in (|1()|) . This allows to 
adopt the same time-stepping method, as will be specified in subsection 12.31 

To quantify the spectral-space effect of multiscale forcing, and also to be 
able to concisely formulate the different forcing procedures in the next subsec- 
tion, we consider the kinetic energy. The equations which govern the Fourier 
coefficients (|10|) can be written in index notation as: 

^ ''^^) " ^ ^^^^ 

where ^o(k, t) = Z)a/3VF/5(k, t) is the nonlinear term. Multiplying this equation 
by the complex-conjugate u*(k, t) and summing over the three coordinate 
directions, we obtain the kinetic energy equation: 

+ 217^^ E{k, t) = <(k, t)^^{k, t) + <(k, t)F,(k, t) (14) 

where Eik, i) = ^ |u(k, t)\^ is the kinetic energy in mode k. Introducing the no- 
tation for the rate of energy transfer T(k, i) = n* (k, t)^Q,(k, i), the rate of 
energy injection by the forcing Tp(k^t) = u* (k, t)FQ,(k, i) and the energy dis- 
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sipation rate e(k, t) 



2i'k'^E{k.,t), we can write H14() equation as: 



dE{k,t) 

Ft 



e{k,t) +T(k,t) +TF{k,t) 



(15) 



This formulation clarifies that the rate of change of kinetic energy -E(k, t) is 
connected with dissipation, expressed by the viscous term e(k, t), with transfer 
to/from different wavenumbers, expressed by T(k, t), and with the forcing term 



The different contributions to the rate of change of the kinetic energy typi- 
cally act in distinct wavenumber regions. The forcing term Tp(k, t) is non-zero 
in the forced modes only. In this paper the collection of forced modes will al- 
ways contain a low wavenumber band corresponding to large-scale forcing of 
the flow. In addition, possible higher wavenumber contributions can be in- 
cluded in Tp(k, t). In contrast, energy dissipation e(k,t) is defined in the en- 
tire spectral space, but it is dynamically important primarily for the high 
wavenumber range, i.e., acting on structures below the dissipation length- 
scale. Finally, the transfer term r(k, t) is basic to the development of an en- 
ergy cascade and is a dominant contribution for wavenumbers in an inertial 
range [33]. In the multiscale forcing cases, we will also introduce forcing gener- 
ally in the same region as where the transfer r(k, t) is dynamically important. 
Hence, the effects of the multiscale forcing relate directly to the 'competition' 
between the dynamics introduced by the forcing procedure and the 'natural' 
transfer of energy to other modes in the spectrum. 

In the formulation of forcing procedures and in the evaluation of the kinetic 
energy dynamics, one frequently adopts shell-averaging. The basic operation 
consists of averaging over spherical shells of thickness 27r/ Lb centered around 
the origin. The n-th spherical shell is given by |^(n — 1/2) < |k| < ^{n + 1/2) 
and will be denoted by ]K.„. Applying shell-averaging to a function /i(k, t) 
defined in spectral space we obtain: 



where P„ is the number of modes in the n-th shell. Applying the shell-averaging 
l|16|) to the energy equation (|15|) we end up with: 



T>(k,i). 




(16) 



dE{n,t) 
Ft 



e{n,t) +T{7i,t) + TF{n,t) 



(17) 



which indicates that the interpretation of the various contributions to the rate 
of change of the kinetic energy at mode k also applies to the shell-averaged 
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formulation. In literature it is common to introduce a numerical correction 
factor when averaging over shells. This is used to compensate for the nonuni- 
form distribution of modes within the discrete spherical shells [11,26]. We will 
follow the convention used in [11,51,52] when presenting the energy-spectra. 
This implies that we multiply h{n, t) by a factor 47rn^ which is associated with 
the 'expected number modes' within the discrete shell. The definition of the 

energy spectrum that we will adopt is given by En = ^47rn^/P„^ -^'(k, t). 

Finally, summing H15|) over all wavevectors k or, equivalently, ()17() over all 
shells yields the evolution equation for the total energy in the system: 

^ = -e{t) + fp{t) ■ h{t) = J2^Pnh{n,t) = Y,^h{k,t) (18) 

where use was made of the fact that the contribution of the transfer term 
r(k, t) is such that it only re-distributes energy over the various modes, which 
implies that its sum over all wavenumbers T{t) = 0. 

Next to spherical shells, it is convenient to introduce spherical wavenum- 
ber bands which consist of several adjacent shells. We denote the wavenum- 
ber band which consists of ^{n^ — 1/2) < |k| < -|- 1/2) by IK^.p, where 
m < p. The corresponding average over IK^.p of a function /i(k, t) is given by: 

1 ^p^^(^^^)^ 1 ^^(k,t) ; = X] P„ (19) 

■^™>P n=m ^""'P K,„,p n=m 

To complete the computational model, we will next introduce the explicit 
forcing strategies that will be investigated in this paper. 



2.2 Explicit forcing procedures 

Forced turbulence in a periodic box is one of the most basic numerically 
simulated turbulent flows. It is achieved by applying large-scale forcing to 
the Navier-Stokes equations. As a result, at sufliciently high Reynolds number 
the well-known turbulent cascade develops in an inertial range of scales which 
are much smaller than the length-scale of the forced modes [27,28,33]. The sta- 
tistical equilibrium that is reached is characterized by a balance between the 
input of energy through the large-scale forcing and the viscous dissipation at 
scales beyond the Kolmogorov dissipation scale. 

Various forcing procedures have been proposed in literature. Generally, 
if the forcing is restricted to large scales only, the specific details of the proce- 
dure do not have such a large effect on the properties of the developing inertial 
range at sufficiently finer scales. However, since we wish to extend the forcing 
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to act on a wide range of scales simultaneously, including parts of an iner- 
tial sub-range, the differences between alternative forcing procedures become 
more pronounced. Investigating these differences is an essential step toward 
quantitative modeling of flow through complex gasket structures and forms 
the main focus of this paper. In this subsection, we will recover the definition 
and some of the motivation for several characteristic forcing procedures. 




Figure 2. Definition of two-band forcing in spectral space (a) and locaUzation of forcing within 

a slab in physical space (b). 

In multiscale forcing, the flow is agitated over a wide range of modes. To in- 
vestigate the effects of such forcing we will focus on cases in which one addi- 
tional spherical band of scales is forced, next to the common forcing of the large 
scales. We consider the general situation as depicted in Fig. Efa). The large 
scales are in the range k < ko and an additional band of small scales is defined 
by fci < k < k2- The forcing method can also incorporate cases in which only 
part of the domain is occupied by a complex obstruction, as sketched for the 
case of a slab in Fig.|2l^b). In fact, by introducing an 'indicator' function 0(x, t) 
to locate the complex object within the flow domain (0 = outside the re- 
gion occupied by the object and 1 elsewhere), the forcing can accommodate 
such spatial localization in a flexible manner. In spectral space, the introduc- 
tion of 0(x, t) implies that the forcing term in spectral space is represented by 
the convolution product of the actual forcing -F(k, t) and the Fourier-transform 
of the indicator function. However, in the present paper such complications 
will not be included and we will only consider forcing procedures applied in 
the entire physical domain. 

Forcing procedures may be classified in different ways. We first distinguish 
forcing schemes which keep the total energy in the collection of forced modes 
identical to its value in the initial condition. This will be referred to as class 'A' 
forcing procedures. Next, we identify forcing schemes which are characterized 
by a constant energy input rate, introduced via the collection of forced modes. 
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This group will be referred to as class 'B'. In either class of schemes, the flows 
develop around a well defined statistically stationary state, but time-dependent 
variations in the total energy and in the energy input rate may occur. 

Apart from a distinction concerning the way energy is introduced into 
the flow, one may classify forcing schemes as 'deterministic' or 'stochastic'. 
Stochastic forcing schemes may introduce an element of uncorrelated random- 
ness, e.g., by restricting the forcing to a random subset of the collection of 
forced modes every time the forcing is invoked. However, the primary ques- 
tion of locality of the modulation of the energy spectrum can be addressed 
more directly using deterministic schemes and in this paper we will restrict 
to these procedures. We next introduce some characteristic forcing schemes in 
either class 'A' or class 'B'. 

Class 'A '; constant- energy forcing. Various methods can be formulated 
which are such that the kinetic energy in the forced modes remains constant. 
The simplest possibility arises by requiring that UQ,(k, t) itself remains constant 
for all k in the collection of forced modes. This was first proposed in [43] and 
implies for the forcing in spectral space: 

Al : ^^(k, t) = uk'^UaiK t) - ^'^(k, t) (20) 

One readily verifies, using (\l'6\i . that dtUa(}^,t) = and in particular this 
implies that 5t£'(k, t) = for each of the forced modes. Hence, also the total 
kinetic energy contained in all the forced modes stays constant in time. The en- 
ergy input rate corresponding to (^0]) is given by Tp(k, t) = e(k, t) — r(k, t) 
for each of the forced modes. This input rate may vary considerably in time, as 
the unsteady flow will lead to a strong time-dependence of the energy transfer 
T for the forced modes. 

The basic method (|2()j) has motivated the formulation of a number of ex- 
tensions which relax the requirement that the Fourier-coefficient is strictly 
constant. In [8] the method was modified to require that |u(k, t)| = const, i.e., 
equal to its initial value, for the forced modes. This allows for the possibility 
that the phases of the Fourier-coefficients may evolve in time. The correspond- 
ing forcing is given by: 

^-C'- *> = (-'^ - «»(^- ') = (^(^- ') - ^(^- '» ilil <^'' 

One may readily verify that this implies 9f£^(k, t) = for the forced modes. 
Forcing expressed in l\'2U\i or (PT|) was found to yield quite large fluctuations 
in the energy input rate [45] . 
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Typically, the forced modes are ordered according to the wavenumber shell to 
which they belong. A shell-oriented simplification of ()21() was proposed [25,26]: 



F^{k,t) = {e{k,t)-Tik,t))^^M. (22) 

This forcing also preserves the total kinetic energy in the forced modes. 
The three forcing procedures (|2()j) . (|^T|) and are quite comparable, both 
in terms of their fluid-physics motivation and in terms of their turbulent flow 
predictions. Therefore, we will only present actual simulation results obtained 
with (|20j) . which are quantitatively representative for the other two forcing 
procedures in this group. 

The forcing methods described so far preserve the kinetic energy that is 
contained in the collection of forced modes. However, considerable variations 
in the total energy in the system can still arise. The reverse can also be real- 
ized, i.e., forced turbulence in which the total kinetic energy in the system is 
constant, but the energy in different modes may vary in time. For this pur- 
pose, the forcing should not be formulated in terms of quantities related to 
individual modes or shell-averaged values, but rather contain averages over all 
modes [18,30]. The case of forcing in a single shell with P modes can readily 
be specified. Specifically, if we replace the shell-average (•) in the amplitude 

factor in H22|) by the average over all modes (•) and use the fact that T = 0, 
we obtain the forcing 

The A2-forcing implies an energy input rate Tp = e{t) and thus by ()18|) 
dE/dt = 0. This method corresponds exactly to the negative viscosity pro- 
cedure used to maintain quasi-steady turbulence direct numerical simulations 
results reported in [21,23,24,50]. Extension of A2-forcing to multiple shells 
can be realized in a number of ways. This will be described in more detail 
momentarily. A2-forcing will be compared to Al-forcing in the next section. 



Class 'B': constant- energy-input-rate forcing. Next to forcing methods 
that can be associated with constant-energy, one may define forcing procedures 
in which the total energy input rate Tp is constant. We first present such 
forcing methods with reference to a single band of forced modes. The way 
in which the energy input is distributed over several bands will be specified 
afterwards. 
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A central example in the class of constant-energy-input-rate forcing methods 
was presented in [16]. Changing e{t) in H23|) into the constant energy input- 
rate e^i; tlie corresponding forcing term may be written as 

The energy input-rate is found to be Tp{t) = e^, as desired by construc- 
tion. The total energy in the system is no longer constant but governed by 
dE{t)/dt = — + £w which implies that the statistically stationary state 
that develops will show a dissipation rate that fluctuates about Syj. This type 
of forcing was also studied in [7,36,52]. Further extensions of the basic forcing 
procedure H24|) can be proposed in which an extra factor k~'^] q > arises in 
the definition of F^. Such an extra factor implies that the forcing of higher 
wavenumber shells can be made to correspond to a specific shape (usually 
A:-5/3 to more directly 'impose' Kolmogorov turbulence). These forcing proce- 
dures will not be considered in this paper; for further details see [9,37]. 

Similar to A-forcing methods, one may formulate related procedures which 
are defined in terms of shell-averaged quantities. For example, analogous to 
(|22j) . we may replace i?(k, t) in (|24j) by E{n,t) to define the forcing of modes 
in the n-th shell. This type of forcing was found to yield basically the same 
results as those based on (|24|) and will not be presented explicitly in the rest 
of this paper. 

The final forcing procedure that we will include in this paper was proposed 
recently in [32] . It was motivated as a model of flow through a fractal gasket 
which functions as a multiscale stirrer. This particular forcing may be asso- 
ciated with a constant energy input rate for the entire system. We modify 
the original forcing procedure slightly and considered in particular 

B2: F„(k,t) = - J e^iKt) (25) 

where K. denotes the set of forced modes. In this formulation, the complexity 
of the object is parameterized by the exponent (3 which is related to the fractal 
dimension Dj of the object through (3 = Dj — 2. The vector e(k, t) has the 
form: 

u(k, t) k X u(k, t) , , 

|u(k,t)| |k||u(k,t)| 



which contains a part in the direction of u and a part that is perpendicular to u. 
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Since eo = |u| = y^2E(k, t) we find for the energy input rate: 

In contrast to Bl-forcing in which the energy input rate is constant in time for 
each of the forced modes separately, this 'fractal forcing' procedure only implies 
a constant energy input rate for the entire system. In fact, after summation 
over all forced modes the total energy input rate is found to be equal to 
Tp(t) = Ew Correspondingly, we find for the evolution of the total kinetic 
energy dE/dt = —e{t) + Ey^, i.e., identical as obtained before for Bl-forcing. 
In the original formulation in [32] the energy input rate was replaced by 
the total dissipation rate which implies that E is constant in time. 

So far, the Bl- and B2-forcing methods were defined with reference to a sin- 
gle band of modes. This band was assumed to contain P modes and was 
identified by K. The total energy input rate £w was available to this band. In 
case more bands are forced simultaneously, the way the energy input-rate is 
divided over the individual bands, and among the modes within each band, 
needs to be specified. For two forced bands and IKma.pa with Pm.i,pi 

and Pm2,P2 modes respectively, such a partitioning involves two steps. First, 
a fraction e^^i = ae^ of the total energy input-rate is 'allocated' to the first 
band and the remainder 6^,2 = (1 — CL)ew is used in the forcing of the second 
band (0 < a < 1). Second, we divide the energy input-rate that is available for 
each band equally over all modes in the corresponding band. As an example, 
two-band Bl-forcing may be defined as 

■ "^'"^■'' = P„.„2i;(k,,) • '^'■> 

(1 -a)g^ Ua(k,t) I. ^ 1.^1. /oox 

= ; otherwise 

The two-band formulation of B2-forcing can be specified analogously, replac- 
ing Ew by either aSw or (1 — a)ew and K. by Kmi,pi or Km2,p2i respectively. 
Extending A2-forcing to more bands can be done in a similar way in which 
a fraction ae{t) is associated with the large-scale band and the remainder with 
the second band. The specific choice of Pmi,pi and Pm2,p2 above implies that 
the energy is equally distributed between all modes within a forced band. We 
can go one step further and require the equal distribution of Ew over the forced 
shells contained in the bands. This implies changing Pmi,pi and Pm2,p2 ™to the 
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number of modes P„ for each forced shell. Extension to more forced bands can 
be formulated analogously. 

In the next subsection the numerical method will be specified in some detail, 
to complete the description of the physical and computational modeling used 
in this paper. 



2.3 Computational method and parallelization 

In this subsection we first specify the time-stepping method, then sketch 
some aspects of the implementation and subsequently describe the validation 
of the method. 



Time evolution. To simulate the spectral solution governed by equations 
()1U() and H12|) we first rewrite these equations in a more general form having 
in mind that the evolution due to the diffusive terms can be computed exactly 
by introducing integrating factors e'^'^ * and e'^^ *, respectively [6]. In fact, (jlOj) 
and ((T^ may be expressed as: 



gU(k,t) 
di 



G(U(k,t)) 



(29) 



where 



U 



u(k, i)e'^'='* 
c(k,f)e'^'='* 



(DW(k,i) +F(k,t))e'^*^'* 
Z(k, t)e'^'='* 



We use a constant time-step At to obtain the solution at times tn = to + nAt. 
A four-stage, second-order, compact-storage Runge-Kutta method was imple- 
mented. The advancement of the solution over a full time-step requires four 
steps of the form: 



U(k, tn+j) = U(k, tn) + 7 At G(U(k, (30) 

The intermediate solutions in the different stages can be found as follows. 
In stage 1 we adopt (7,^) = (1/4,0), stage 2 requires (7,^ = (1/3,1/4), 
stage 3 uses (7,^) = (1/2,1/3) and stage 4 completes the step with (7,^) = 
(1,1/2) [15]. 

We consider turbulence in a cubic box of side Lj, with periodic boundary 
conditions and assume that the flow is statistically isotropic which implies 
that we require the same resolution in each coordinate direction. The direct 
numerical simulations will employ a resolution of N^, where is the num- 
ber of spectral-space grid-points that is used in each direction. This restricts 
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the set of wavenumbers to tIq, = 0, ±1, ±2, . . . , it — 1) , —N/2. The cut- 
off wavenumber is given by /cmax = vrA^/L;,. In physical space this corresponds 
to a uniform grid = jL^/N, where j = 0,1,2, N — 1 in each coordinate 
direction. 

We use the pseudo-spectral discretization method, i.e., the spatial derivative 
terms in the Navier-Stokes and passive scalar equations are computed via sim- 
ple multiplications in the spectral space. The nonlinear terms in the equations 
are evaluated in physical space to avoid the evaluation of several computation- 
ally intensive convolution sums [6]. This procedure requires three steps. First, 
the Fourier-coefficients u(k, t) and c(k, t) are used to obtain the velocity and 
scalar fields in the physical space. Subsequently, the velocity-velocity products 
and the velocity-scalar products are determined in physical space and finally 
the associated Fourier-coefficients of these products are obtained. 

The finite resolution may give rise to well-known aliasing errors. In fact, 
the product of two Fourier-series based on a resolution with N points gives 
rise to more small-scale modes than can be supported by the grid. As a result, 
these contributions can appear on the A^-point resolution as seemingly lower 
wave-number modes. A detailed discussion of techniques allowing the partial or 
full removal of the aliasing error can be found in [6]. Several of these methods 
were tested, as described in the appendix, closely following [40,41]. The method 
of two shifted grids and spherical truncation was used in actual simulations. 
This removes the aliasing error completely which was found to be essential, 
especially to maintain the characteristics of the small turbulent scales. 



Data decomposition and fast Fourier transforms. The simulation software 
was implemented in Fortran 90 and parallelized based on the framework given 
in [52] using the Message Passing Interface (MPI) [38]. Data are stored using 
the Hierarchical Data Format (HDF5) [19] which is a file format and library 
designed for scientific data-storage and handling. The choice of HDF5 was 
motivated by the flexible data-exchange between different platforms and its 
support of parallel I/O. High performance computations were done at SARA 
Computing and Networking Services (Amsterdam) on Silicon Graphics (SGI) 
Altix 3700 and Origin 3800 CC-NUMA systems (for more details see [42]). 

The critical performance factors in the parallel implementation of 
the pseudo-spectral discretization method are the domain decomposition 
and the algorithm for the three-dimensional Fast Fourier Transform (FFT). 
These two implementation decisions are essential since they determine almost 
all aspects of the data-exchange between domains and most of the floating 
point operations. It is important to obtain a data decomposition which permits 
for fast transfer of data between processors. To obtain parallel Fourier trans- 
forms we adopted procedures from two libraries: SCSL [17] and FFTW [13]. 
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Moreover, since access to memory and the number and speed of available 
CPU-s may differ considerably among different computational platforms, sig- 
nificant improvements in the processing time can be achieved by platform- 
dependent optimization. 

The speedup of the parallel implementation was checked by simulating de- 
caying turbulence at a resolution of 256'^. A time-interval < t < 0.05 was 
considered. This case corresponds to 28 time-steps with 5 data evaluation and 
reporting stages. In a non-dedicated SGI Altix 3700 environment we obtained 
on 4, 8, 16, 32, 64 processors the following speedup numbers: 3.9, 7.5, 14, 26, 47, 
respectively. The best performance results were obtained by a cache-unfriendly 
parallelization along the second array dimension. This gives the opportunity 
of minimal data exchange and reshuffling between processors and illustrates 
that the speed of the processors overwhelms the abilities of direct access to 
the memory. This was found to be the critical issue for the hardware that was 
available. 

Code validation. To validate the implementation of the pseudo-spectral 
method, decaying homogeneous isotropic turbulence was simulated at two dif- 
ferent Reynolds numbers. The initial condition was taken from [35], which 
was generated on the basis of the Pao spectrum [39] . For further details we re- 
fer to [34]. This flow was studied extensively using high-order finite- volume 
discretization and explicit Runge-Kutta time-stepping. Special attention was 
given to the degree of convergence that could be achieved using the finite 
volume approach. These data provide a clear point of reference with which 
the present pseudo-spectral flow-solver can be compared. 
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Tabic 1. The value of kmaxV assoeiated with different resolutions. The Kolmogorov scales are 
ri = 5.87 • 10"^ and tj = 2.07 ■ 10~^ for Rx = 50 and Rx = 100, respeetively. 



A first, global assessment of the resolution that is achieved may be inferred 
by evaluating the product of the cut-off wavenumber and the observed Kol- 
mogorov dissipation length-scale r] = L(3i?|/20)^'^/^ in terms of the Taylor 
Reynolds number Rx computed for the initial condition (see H33() for the defini- 
tion) and integral length L = 1/2. In order to resolve all dynamically relevant 
length-scales, including the dissipation length-scale it is required that /cmax?? 
is sufficiently large. A commonly accepted criterion of adequate spatial reso- 
lution is that /cjnax?? > 1- When the focus is on higher-order statistics, it is 
preferred to use larger values {k^anxV > 3/2) [11,46]. In Tabled the values of 
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kma^V are presented for the two computational Reynolds numbers considered 
Re = 1060.7 and Re = 4242.6 which correspond to Rx = 50 and Rx = 100. 
We observe that in the first case a resolution of at least 64"^ is required to 
achieve full resolution, while in the second case the minimal required resolu- 
tion moves up to 192'^. 

For validation of the code, the flow was simulated for more than two eddy- 
turnover times and a number of quantities were monitored: 



Rx{t) = X{t)u{t)/u ; u{t) = 

{{dv,{^,t)/dx^f) 

Jlirj — 57^7 

((9t;i(x, 0/5x1)2)^/2 

The operator (•) in H34|) refers to volume averaging. 



Total energy 
Taylor microscale 

Taylor Reynolds 

Longitudinal skewness 



(31) 
(32) 



-E{t) (33) 



(34) 




Figure 3. Prediction of total energy E, Taylor microscale A, Taylor Reynolds number Rx and 
longitudinal skewness Si at an initial = 50 (a) and Rx = 100 (b) with a 
finite-volume [35] (solid) and the present pseudo-spectral (dotted) code. 

In Fig. 131 a comparison is made between simulation results obtained with 
the pseudo-spectral method at = 512, and with the high-order finite- volume 
discretization method [35]. For each of the quantities an almost perfect agree- 
ment may be observed. In Fig.^we assessed the convergence of the predictions 
as function of the spatial resolution. In this figure we replaced the longitudinal 
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skewness Si by the skewness: 



For homogeneous isotropic turbulence the value of S should be equal 
to 0.5 [2], which is quite well approximated in the simulations. This quantity is 
quite sensitive to the spatial resolution and is therefore a good indicator of ap- 
propriate spatial resolution. We observe that the different predictions display 
a clearly distinguishable convergence toward the grid-independent solution. 
Specifically, results obtained for resolutions higher than 64^ at R\ = 50 and 
192^ at Rx = 100 are quite indistinguishable, consistent with the criterion 

that /Cmaxf? > 1- 

In the next section we turn to the effects that different multiscale forcing 
procedures have on the developing turbulent flow. We will focus in particular 
on the modifications that arise in the kinetic energy spectrum. 



3 Modulated cascading by broad-band forcing 

The explicit forcing in different wavenumber bands can have a strong effect 
on the developing turbulent flow. We discuss the modifications of the energy 
spectrum arising from 'constant-energy' (class 'A') or 'constant-energy- input- 
rate' (class 'B') procedures. The various forcing strategies will be shown to 
qualitatively correspond to each other, provided the total dissipation-rate 
and the spectral energy distribution are commensurate for the different class 
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'A' and 'B' forcing strategies. We will specify this inter-relation in more detail 
momentarily. As point of reference, we will first turn our attention to forcing 
of the large scales only. Subsequently, we consider two-band forcing and inves- 
tigate in particular the effects of variation of the strength and location of the 
small-scales band on the developing flow. 

In the sequel, we consider time-averaged properties of the developing turbu- 
lent flow defined by: 



where T is sufficiently large. In all cases to = 5 in order to allow the averaging- 
process to start from a properly developed quasi-stationary state. The aver- 
aging is continued up to T = 25, which corresponds to approximately 40 
eddy-turnover times. This was found to provide an accurate representation of 
the long-time averages, leading to relative errors below 5 percent, measured in 
terms of the ratio of the standard deviation and the mean signal. This proce- 
dure was applied to obtain the time-averaged kinetic energy spectra as well, 
which are very effective for monitoring changes in the kinetic energy dynamics 
due to the forcing. 

Large-scale forcing. To create a point of reference, we first consider forced 
turbulence in which energy is introduced to the system only in the first 
shell Ki^i. We adopt ko = Sir referring to Fig. [21^a) and force all 18 modes in- 
side this band. The computational Reynolds number Re = 1060.7 and the size 
of the computational domain Lf, = 1. The spatial resolution was taken to 
be 128^, which provides ample resolution of these cases, similar to what was 
established in subsection 12.31 

In order to be able to quantitatively compare results obtained with the 
different forcing strategies, care should be taken of properly 'assigning' a level 
for the energy dissipation-rate and the spectral energy distribution. For this 
purpose, we may consider simulations with the A2-method to be central in the 
sense that the other three forcing strategies may be specified with reference to 
it. In fact, if we generate an initial condition with a certain total kinetic energy, 
then A2-forcing yields an evolving fiow which becomes statistically stationary 
after some time, while maintaining the same level of total energy. The A2- 
forced simulation can be used to specify the 'corresponding' class-B forcing 
strategies. In fact, the constant dissipation-rate in class 'B' forcing is taken 
equal to the time-average value of the dissipation-rate that is found from the 
A2-forced simulation, i.e., we adopt = (?)t. This procedure was adhered to 



t 



T 




(36) 



to 
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in all cases presented in this section. Finally, in the developed stages of either 
these A2- or B-forced flows, any instantaneous solution may be used to arrive 
at a full specification of the 'corresponding' Al-forcing. The actual choice of 
this instantaneous solution is arbitrary. However, when comparing simulations 
based on Al-forcing that adopt different realizations of the turbulent flow-field, 
we observed that the statistical properties of all these Al-forced cases were 
the same. 




Figure 5. The evolution of the total kinetic energy E (a) and energy-dissipation-rate e (b) for 
the large-scale forcing: Al (dashed), A2 (dotted), Bl (dash-dotted), B2 (solid). 



The evolution of the total kinetic energy E{t) and energy dissipation rate e{t) 
is shown in Fig.|Hl As initial condition for the A2- and B-forced simulations, we 
adopted the velocity field obtained aX, t = 0.5 from the decaying homogeneous 
turbulence simulation discussed in subsection 12.31 To be able to qualitatively 
compare with the Al-forced flow at a similar energetic level we took as initial 
condition the solution from Bl-forcing at t = 5. The total kinetic energy is seen 
to fluctuate around its long-time mean value (of course, apart from A2-forcing). 
As can be seen, the system rapidly develops into a statistically stationary 
state characterized by the input of energy, its transfer to smaller scales and 
dissipation in the viscous range. In Al-forcing the Fourier-coefficients in the 
forced band are all kept constant, i.e., equal to their initial values. The energy 
in the system fluctuates very significantly, which was considered a disadvantage 
of this forcing in [11]. The energy and dissipation levels in Al-forcing differ 
considerably from those obtained with the other forcing strategies. To compare 
A2-forcing with B-forcing, the energy dissipation rate was taken as £u] = (e)^ — 
0.2. The total kinetic energy for Bl-forcing is seen to fluctuate around the 
constant value associated with A2-forcing. A similar impression is observed 
when use is made of the fractal B2-forcing in which the fractal dimension of 
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the stirrer was taken equal to Df = 2.6 [32] which corresponds to an exponent 
P = 3/5 in (EHJ). 

In general, when applied to the largest scales only, all forcing procedures 
mentioned in subsection l2.2l vield similar results. As a further example, the tails 
of the time-averaged spectra were found to be virtually identical to each other, 
which indicates that the properties of the smaller turbulent length-scales are 
not very strongly dependent on the details of the specific forcing. This was also 
established by various other quantities that were investigated. Specifically, the 
Taylor-Reynolds number Rx for the simulated cases was seen to fluctuate in 
the range between ~ 50 up to ~ 60 for all methods. The time-averaged value 
of the skewness was also investigated and found to be very close to 0.5. This 
indicates that a well developed isotropic flow was attained [2]. 



Two-band forcing. In the simulations that adopt two-band forcing we con- 
sider situations in which we introduce energy into the system in a band con- 
sisting of four shells, next to the already described large-scale forcing in the 
first shell. We first compare the different class 'A' and 'B' forcing strategies, 
within this two-band setting. As second, forced band we consider ]ICi7^20- This 
band corresponds to ki = SStt and k2 = 4l7r in Fig. I2ta) and contains in 
total 17284 different modes that are all explicitly forced. The comparison of 
the different forcing strategies shows that the flow-predictions are qualitatively 
comparable. Subsequently, we therefore focus on the B2-forcing strategy and 
investigate the effects arising from changes in the strength or the location of 
the second forced band. 

Forcing of a second band implies that we need to additionally specify how 
the energy input is distributed over the bands, the shells within the bands 
and, finally, the modes within the shells. The specification of the A2-forcing 
requires the fraction of the energy-input that is allocated to the different bands. 
We consider the case in which a = 1/5 in (|28j) which corresponds to equi- 
partitioning of the energy- input over the five shells that are forced. The forcing 
within the second band is further specified by assigning an equal energy-input 
rate to each of the four shells contained in it. Finally, each of the modes 
in a particular shell n receives an equal share of the energy-input to that 
shell, taking the number P„ of modes in the particular shell into account. To 
compare the 'A' with 'B' forcing strategies we adopt the same method as above 
for specifying the parameter e^. Specifically, the total energy injection for the 
'B' methods was given as = = 1, in terms of the time-average of the 
total dissipation rate in the A2-forcing. Moreover, the same equi-partitioning 
of the energy-input as in A2-forcing was adopted. Finally, the Al-forcing is 
derived from the field that was obtained at t = 5 with the Bl-forcing. We 
verified that the statistical properties of the Al-forced flow are insensitive to 
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the particular choice of initial field used to define this forcing method. 



0.5 1 • • • ■ 1 :' 




Figure 6. The evolution of the total kinetic energy E (a) and cncrgy-dissipation-rate e (b) for 
two-band forcing: Al (dashed), A2 (dotted), Bl (dash-dotted), B2 (sohd). 



As may be noticed by comparing Fig. [S] with Fig. 1^1 the two-band forcing 
leads to a strong increase in the total energy dissipation-rate, while the to- 
tal kinetic energy present in the flow is quite unaffected by the second forced 
band. The increase in the dissipation-rate is particularly strong for Al-forcing. 
Hence, the high-A: forcing changes mainly the distribution of energy over 
the scales and not so much the actual energy content. By changing the strength 
and location of the forcing, we have the possibility to control, and to some 
extent manipulate the way the energy is distributed and hence indirectly in- 
fluence the large- and small-scale transport properties of the flow. We turn to 
this aspect next, by focusing explicitly on the kinetic energy spectrum. 

The compensated kinetic energy spectra E{k) = (e)^ "^^^ k^/'^ {Ek)t that are 
obtained with the different two-band forcing methods are collected in Fig. EJ 
The modiflcations in the spectrum, relative to the case of large-scale forc- 
ing only, are localized primarily in the region close to the forced band. All 
forcing methods are seen to yield qualitatively quite similar results. Next to 
the expected modifications near the explicitly forced band, we observe that 
the two-band forcing also affects a much wider set of larger scale modes. In 
fact, a significant depletion of the kinetic energy in a range of scales 'ahead 
of the forced region, is readily appreciated. This indicates that the agitation 
of a small band of modes can induce large changes in a rather wide part of 
the spectrum which further characterizes the type of turbulence-control that 
one may achieve with explicit forcing. 

We next turn to the second part of this section and consider the effects 
of varying the spectral support and the strength of the second forced band. 
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Figure 7. Compensated energy spectrum for two-band forced turbulence (fc < Stt and 
337r < fc < 4l7r) with different methods: Al (dashed), A2 (dotted), Bl (dash-dotted), B2 (soUd). 




Figure 8. Compensated energy spectrum for two-band forced turbulence with the B2-method and 
the same energy inputs em,i = £m,2 = 0.15 to the fc < Stt band and various locations of the second 
band: 97r < fc < 17ir, 17tt '< k < 23-k, SSvr < fc < 4l7r, 497r < fc < STtt (dashed, dash-dotted, o, o). 
The spectrum obtained with large-scale forcing at e-w = 0.15 in fc < Stt band (solid). 



The qualitative similarity of the different two-band forcing methods as seen 
in Fig. [3 allows to concentrate on only one of the forcing methods. We adopt 
B2-forcing in the sequel. In Fig. |H1 we illustrate the effect of variation of the 
spectral support of the second band. Relative to the case of large-scale forcing 
only, we observe that the tails of the spectra are quite unaffected. However, 
the injection of energy in the second band is seen not only to increase the en- 
ergy in the forced scales but also to deplete the energy in all the larger scales. 
Moreover, the 'up-scale' effect of energy-depletion is more pronounced in case 
the second band is moved toward smaller scales. 
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Figure 9. (a) Compensated energy spectrum for two-band forced turbulence with the B2-forcing 
method and energy input rate e^^i = 0.15 in fc < 3tt band for different strengths of forcing in the 
337r < fc < 4l7r band: e„,2 = 0.075,0.15,0.30,0.45,0.60,0.75,0.90 {□, dotted, dashed, dash-dotted, 
>, o, o). Large-scale forcing with ew = 0.15 in the fc < 37r band is denoted by the solid line, (b) 
Corresponding time-averaged total kinetic energy with standard deviations. 



The control over the flow that is available with two-band forcing is examined 
further by investigating the effects of varying the strength of the forcing in the 
second band. We kept the energy input rate for the first A; < Svr band equal 
to £w = 0.15 and varied the strength of forcing in the second SSvr < k < 4l7r 
band adopting = 0.075 . . . 0.90. The corresponding compensated energy 
spectra from these simulations are shown in Fig.|^a). We observe that a higher 
energy input-rate in the second band leads to a more pronounced peak in the 
spectrum which shifts to lower values of ki] with increasing of the second 
band. Simultaneously, the value of E{k) decreases for the low-A; modes with 
increasing e^. This forcing of the second band allows to quite independently 
control the spectrum, at roughly the same total energy content in the flow. 
In fact, variation of of the second band by a factor of about 10 is seen to 
lead to a comparably strong increase in the peak value of the spectrum while 
the total energy level {E)t is increased by only 15% as seen in Fig. Efb). 

In the next section we will examine how the changes in the flow properties 
due to the two-band forcing in spectral space influence the physical space mixing 
efficiency of a passive scalar. 



4 Small and large scale mixing efficiency 

The consequences of explicit broad band forcing not only express themselves 
in modulated energy cascades. The mixing properties of the evolving turbulent 
flow in physical space also depend significantly on the forcing that is applied. 
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In this section we quantify the mixing-efficiency by monitoring geometric prop- 
erties of evolving level-sets of an embedded passive scalar. The numerical inte- 
gration method that is used to determine these level-set properties is described 
in subsection 14.11 The ensemble-averaged simulation results are discussed in 
subsection 14.21 we establish to what extent the two-band forcing can be used 
to control the maximal rate of mixing and the total accumulated degree of 
mixing. 



4.1 Level-set evaluation to quantify mixing 

To illustrate and quantify the influence of two-band forcing on the turbulent 
dispersion of a passive scalar field we analyze the evolution of basic geomet- 
ric properties of its level-sets. As a result of the turbulent flow these level- 
sets become highly distorted and dispersed across the flow-domain. Specifi- 
cally, we concentrate on the surface-area and the wrinkling of these level-sets. 
We adopt a specialized integration method to determine these geometric prop- 
erties, as developed in [14]. This method is based on the Laplace transform 
and avoids the explicit construction and integration over the complex and 
possibly fragmented scalar level-sets. With this method an accurate and effi- 
cient evaluation of the evolving mixing efficiency can be achieved which allows 
to quantify the increased complexity of the flow in relation to the two-band 
forcing that is used. 

Basic geometric properties of a level-set S{a,t) = {x G M | C(x, t) = a} of 
the scalar C{x,t) may be evaluated by integrating a corresponding 'density 
function' g over this set. In fact, we have: 

Ig{a,t)= [ dAg{x,t)= [ dx5{C{x,t)-a)\VC{x,t)\g{x,t) (37) 

Jsia,t) Jv 

where the volume V encloses the level-set S{a,t) [31]. Setting g{x,t) = 1, 
g{:>c,t) = V-n(x, t) or g{x,t) = |V-n(x, t)|, we can determine the global 
surface-area, curvature or wrinkling of S{a, t) respectively. Here 
n(x, t) = VC(x, f)/|VC(x, f)| is a unit normal vector, locally perpendic- 
ular to the level-set. The divergence of this vector-field is directly related to 
the local curvature of the level-set. 

We will focus on the evolution of the surface-area A and the wrinkling W. 
The scalar C is scaled to be between and 1; we will primarily consider 
the level-set a = 1/4. In particular we monitor: 
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By determining -dA and 'dw we may quantify the rate at which surface-area 
and wrinkUng develop, the maximal values that are obtained and the time- 
scale at which these are achieved. The corresponding cumulative effects are 
given by 



These cumulative measures express the total surface-area and wrinkling that 
has developed in the course of time. In particular applications, e.g., involving 
combustion in diffusion flames, the cumulative surface- area and wrinkling ex- 
press the total 'chemical processing capacity'. Here, we will determine these 
cumulative effects in order to characterize the different two-band forcing pro- 
cedures. 

To establish the influence of forcing on turbulent mixing properties we sim- 
ulated the spreading of a passive tracer at Schmidt number Sc = 0.7. The sim- 
ulations were started from a spherical tracer distribution of radius r = 3/16. 
The scalar concentration was set equal to 1 inside this sphere and outside. 
A localized Gaussian smoothing of this C-distribution was applied near the 
edge of the initial sphere to avoid resolution problems. The fractal forcing B2- 
procedure as defined in (|25|) was adopted. We performed numerical simulations 
in which the energy input-rate and the spectral support of this two-band 
forcing were varied. 

As point of reference we adopted large-scale forcing in the ICi^i shell with 
an energy injection-rate = 0.6. The resolution requirements were satisfac- 
torily fulfilled: k^axV ranges from 2.3 to 3.5 using a resolution in the range 
128^ - 192^ grid-cells. For the passive scalar these resolutions correspond 
to A:jnax'?oc in the range from 3 to 4.5, where r]oc is the Obukhov-Corrsin 
scale [47]. To study the influence of two-band B2-forcing we applied supple- 
mentary forcing either in a region situated near the largest scales of the flow, 
i.e., IKs^g or further separated, i.e., Kia^ig. In case two bands are forced, the en- 
ergy input-rate for the IKi^i shell is e^j^i = 0.15, while the second band is forced 
using ew,2 = 0.45. In this way the total energy level is kept at comparable lev- 
els in the different cases. A qualitative impression of the effect of these forcing 
procedures may be observed from the snapshots shown in Fig. 1101 The veloc- 
ity and passive scalar display considerably more small-scale features in case of 
two-band forcing, particularly in case of high-/c forcing. To quantify this qual- 
itative impression we apply the level-set analysis discussed above. The results 
will be presented in the next subsection. 




(39) 



28 



Mixing in manipulated turbulence 




Figure 10. Snapshot of vertical velocity field iso-surfaces (left) and passive scalar concentration 
(right) at t = 0.5 for large-scale forcing Ki,i with = 0.6 (a), or with euj,i = 0.15 in the first shell 
and complementary forcing ew,2 = 0.45 in Ks^g (b) or Ki3,i6 (c). In the velocity field snapshots the 

red iso-surface corresponds to U2 = 0.2 and the blue iso-surfaces to U2 = —0.2. The iso-surface at 
C = 0.25 is shown for the passive scalar. 

4.2 Surface- area and wrinkling 

In this subsection, we compare instantaneous and accumulated mixing prop- 
erties for large-scale forcing and different two-band forcing. The total energy 
input rate to the flow is kept constant at 0.6; a fraction e^^i is allocated to the 
first shell and 6^,2 to the second band such that e^,! + £w,2 = 0.6 and e^,! is 
varied from 0.05 up to 0.6. The characterization of the mixing-efficiency was 
based on averaging 20 simulations, each starting from an independent real- 
ization of the initial velocity field. The different initial conditions were each 



Mixing in manipulated turbulence 



29 




Figure 11. Evolution of decaying passive scalar growth parameters: a) surface-area i?^, 
b) wrinkling 'd\Y, c) accumulated surface-area (^^, d) accumulated wrinkling Large-scale 
forcing Ki i with Sw = 0.60 (solid) and complementary two-band forcing (e^,! = 0.15 and 
£w,2 = 0.45) in K5,8 (dashed), Kis^ig (dash-dotted). 



separated by two eddy-turnover times. 

The instantaneous and cumulative effects arising from both the large-scale 
and the two-band forcing are shown in Fig. 1111 The development of the in- 
stantaneous surface-area and wrinkling is qualitatively similar in each case. 
The concentrated initial tracer distribution is in the first stages primarily dis- 
persed by convective sweeping in the turbulent flow. As a result, the level-set 
corresponding to a = 1/4 becomes distorted and both 'Oa and -dw show a rapid 
increase. However, since no source of scalar was included in the computational 
model, molecular diffusion dominates the long-time behavior and leads to "Oa 
and -dw to decrease to zero as t — s- cxd. In between, '&a and 'dw reach their 
maximum. The rapid initial growth is also clearly expressed in Fig. Illf c-d). 
In addition, the cumulative measures and CvK show a clear saturation as 
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Figure 12. (a) Evolution of the decaying passive scalar accumulated surface-area parameter 
for two-band forcing (Ki i - Kis.ie) with different equi-partitions of energy between bands: 
(0.60 - 0.00) (o), (0.45 0.15) (solid), (0.30 - 0.30) (dashed), (0.15 - 0.45) (dot-dashed), 
(0.05 - 0.55) (dotted), (b) Total surface-area and wrinkling at t = 2 for different ew,2 in 
Ki3,i6 (results normalized by the total surface-area and wrinkling for the large-scale forcing). 

We observe from Fig- llir a-b) that forcing of the large scales only, creates the 
highest growth-rate of surface-area and wrinkling. The surface-area reaches its 
maximum value both sooner and at a higher value in this case. In the initial 
stages convective spreading of the tracer dominates over the decay caused by 
molecular diffusion; hence in these stages the agitation of the larger scales plays 
a crucial role in the evolution of the surface-area. The higher band forcing needs 
to compete more directly with the viscous effects and does not induce very 
strong sweeping motions over large distances. Correspondingly, high-A; forcing 
is found less effective in producing surface-area. The more localized distortions 
of the scalar level-sets, as expressed by the development of the wrinkling, are 
less affected by the competition with viscosity, as seen in Fig. Illf b). 

The interpretation of the effectiveness of the mixing in relation to the spe- 
cific forcing alters if we compare the accumulated values for surface-area and 
wrinkling. As may be appreciated from Fig. llir c-d). a significant enhancement 
of the accumulated long-time surface-area and wrinkling arises as a result of 
the explicit agitation of the smaller scales in the flow. Evidently, forcing of 
the smaller scales does not yield a more intense mixing, judging from the in- 
stantaneous values, but does yield an increase in the total surface-area and 
wrinkling, accumulated over time. 

To measure the influence of variations in the strength of the forcing in 
the high-A; band, we focus on the Ki^i and IKia^ig two-band forcing. In par- 
ticular, we keep ew,i + £w,2 = 0.6 and vary the values of 8^,2- The effects on 
the cumulative mixing-efficiency are shown in Fig. I12f a). We observe that an 
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increase in e^^2 implies a slight decrease in the initial growth rate of instan- 
taneous surface-area, but an increase in the long-time cumulative effect. The 
dependence of the long-time cumulative effect on £^^2 is clarified in Fig. ll2f b). 
These simulation results establish the degree of control that may be achieved 
with two-band forcing and the feasibility of such computational modeling. This 
approach may help to identify optimal stirring procedures to which future re- 
search will be devoted. 



5 Summary and concluding remarks 

Various deterministic forcing methods that perturb a turbulent flow in a cho- 
sen range of length-scales were examined. The presented modeling framework 
incorporates the explicit forcing as an integral part. We have shown that with 
a relatively simple forcing model basic properties of complex flows can be 
captured. For example, an enhancement of the energy dissipation by small- 
scale forcing was seen to produce so-called spectral short-cut features, quite 
similar to what was observed experimentally in flow over canopies [12] where 
the kinetic energy is immediately transferred to the smallest scales of the flow. 

Forcing methods agitating the flow in a wide range of scales induce signifi- 
cant differences in the developing flow, compared to the case obtained classi- 
cally in which only the large scales are forced. Various forcing methods were in- 
troduced and shown to produce qualitatively similar results, provided the forc- 
ing parameters correspond to turbulence at comparable total kinetic energies. 
We classified the methods according to constant-energy or constant-energy- 
input-rate and examined these procedures by simulating forced turbulence 
with energy injected in two different bands. The modulation of the turbulent 
flow was investigated for various locations of the second high-fc band in spec- 
tral space. It was shown that the forcing in the second band induces nonlocal 
modulation of the energy spectrum. This was further examined by simulations 
done with different strength of forcing in the high-Zc band controlled by the 
energy injection rate. 

We devoted special attention to a recently proposed multiscale forcing that 
models a flow under the influence of an additional perturbation by a multiscale 
object [32]. We performed numerical simulations of the dispersion of a passive 
scalar field in a turbulent flow that is driven by such forcing. A level-set integra- 
tion method was adopted to quantify general characteristics of mixing quality 
and efficiency. It was found that broad-band forcing causes additional produc- 
tion of smaller scales in the flow which is directly responsible for the localized 
enhancement of the wrinkling of the level-set. In contrast, the surface- area of 
a level-set of the tracer is found to be mainly governed by convective sweeping 
by the larger scales in the flow and hence it is governed to a greater extent by 
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the energy input-rate allocated to the small-A; range. Future study will include 
the spatial localization of the forcing. This can help to model flows that are 
more closely related to realistic physical situation observed in experiments and 
applications. 
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Appendix A: 

To eradicate the aliasing error we study in more detail (i) the random phase 
shifts method and (ii) the method employing two-shifted grids with spheri- 
cal truncation [6,40,41]. In the first case, the aliasing error is only partially 
removed. With additional truncation of the Fourier velocity field coefficients 
the remaining error can be reduced to O(Af^). In the method employing two- 
shifted grids and spherical truncation, the aliasing error can be fully removed 
from the simulations. This specific approach doubles the computational costs 
and memory requirements, compared to the random shifts method. The well- 
known 3/2 method can be used as well, by going to higher resolution and 
truncating the field. This can be done with the lowest number of operations, 
but has higher memory requirements. 

The aliasing error for higher-resolution runs affects mainly the small-scale 
statistics. This is visualized in Fig. lAll where we have shown the Taylor- 
Reynolds number and the longitudinal skewness for decaying turbulence sim- 
ulations with initial Rx = 100 and two resolutions 128'^ and 192'^. The partial 
dealiasing removes the main aliasing error and with the additional truncation 
reduces it to the accuracy associated with the adopted Runge-Kutta scheme. 
There is a small difference between the full and partial removal of the aliasing 
error for the lower resolution of 128'^, but this largely vanishes for the well- 
resolved 192'^ case. 
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Figure Al. Influence of aliasing error for resolution (S'-b) and 192^ (c-d) on Taylor-Reynolds 

number R\{t) and longitudinal skewness Si{t) at an initial = 100 case. Simulations with 
aliasing error (dotted), partial dealiasing without truncation (dashed), partial dealiasing with 
the truncation (dash-dotted), full dealiasing by two grid shifts (solid). Results for the partial 
dealiasing with truncation (dash-dotted) are almost identical to fully dealiased results (solid). 
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